#include "rk4.h"
#include "func.h"
#include "zadacha.h"

void rk (double t, double th, double* x, std::valarray<double> w, input& data) {
    unsigned m = data.q;
    std::valarray<double> k1(m),k2(m),k3(m),k4(m),dydx(m),xd(m),rez(m);

    std::valarray<double> u(data.num_upr);

    for (unsigned i=0; i<m; i++) xd[i] = x[i];
    u = upr(t,w,data);
    dydx = right_parts(xd,t,u,data);
    for (unsigned i=0; i<m; i++) k1[i] = th*dydx[i];


    for (unsigned i=0; i<m; i++) xd[i] = x[i]+k1[i]/2;
    u = upr(t+th/2,w,data);
    dydx = right_parts(xd,t+th/2,u,data);
    for (unsigned i=0; i<m; i++) k2[i] = th*dydx[i];


    for (unsigned i=0; i<m; i++) xd[i] = x[i]+k2[i]/2;
    u = upr(t+th/2,w,data);
    dydx = right_parts(xd,t+th/2,u,data);
    for (unsigned i=0; i<m; i++) k3[i] = th*dydx[i];


    for (unsigned i=0; i<m; i++) xd[i] = x[i]+k3[i];
    u = upr(t+th,w,data);
    dydx = right_parts(xd,t+th,u,data);
    for (unsigned i=0; i<m; i++) k4[i] = th*dydx[i];

    for (unsigned i=0; i<m; i++) {
        x[i] += k1[i]/6 + k2[i]/3 + k3[i]/3 + k4[i]/6;
    }

}
